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Abstract 

Rhythms in electrical activity in the membrane of cells in the suprachiasmatic 
nucleus (SCN) are crucial for the function of the circadian timing system, 
which is characterized by the expression of the so-called clock genes. Intra¬ 
cellular Ca 2+ ions seem to connect, at least in part, the electrical activity of 
SCN neurons with the expression of clock genes. In this paper, we introduce 
a simple mathematical model describing the linking of membrane activity to 
the transcription of one gene by means of a feedback mechanism based on 
the dynamics of intracellular calcium ions. 
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1. Introduction 

Circadian rhythms arise from the cooperative action of a number of en¬ 
dogenous biological oscillators generating daily patterns of many physiolog¬ 
ical and behavioral processes that persist even in the absence of the forcing 
provided by the external light-dark cycle. The master circadian clock in 
mammalians is located in the suprachiasmatic nucleus (SCN), a neuronal 
structure located in the anterior hypothalamus. The cells of the SCN be¬ 
have as a set of cooperative autonomous oscillators with slightly distributed 
frequencies yielding a global circadian (~24-hour period) rhythm. Thus the 
basic oscillatory mechanism leading to the emergence of circadian rhythms 
has an intracellular origin that relies on the negative self-regulation of gene 
expression through transcriptional/translational feedback loops. In the last 
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few years, a number of genes involved in such a regulatory mechanism has 
been identified [I]. 

On the other hand, the neurons in the SON show a characteristic firing 
pattern that changes dramatically with the circadian cycle. They are thought 
to encode the time of the day by adjusting their firing frequency to high rates 
during the day and lower ones at night. Different studies carried out in recent 
times both in Drosophila and in mammals suggest that the electrical activity 
of SON cells provides the driving of the molecular clockwork. Also, keeping 
the SON cells within an appropriate voltage range may be required for the 
generation of circadian rhythmicity of clock gene expression at the single 
cell level [T|. Thus, a fundamental question in circadian biology is how the 
electrical activity may regulate clock gene expression and, conversely, how 
this latter process may alter the electrical activity in SCN neurons. 

The evidence accumulated in recent years suggests that the effect of the 
electrical activity on clock gene expression by SCN cells is probably mediated, 
at least in part, by Ca 2+ ions. In fact, a close relationship between electrical 
activity and Ca 2+ levels has been observed. Resting levels of Ca 2+ in SCN 
neurons exhibit a circadian rhythm that has been detected by using a calcium 
sensitive dye. During the peak in firing at midday, SCN neurons show resting 
Ca 2+ levels of around 150 mM, but these levels drop to about 75 mM during 
times of inactivity. The action potential itself is an important source of Ca 2+ 
in the SCN, regulating Ca 2+ influx into the soma through the opening of 
voltage-sensitive calcium channels [1]. This feature has been shown most 
clearly by a recent work in which a Ca 2+ levels and the firing of the SCN 
neurons were measured simultaneously [2]J. Data from this study show that 
driving the frequency of action potentials in the SCN neurons to 5-10 Hz 
(daytime levels) induces a rise in somatic Ca 2+ levels. This effect can be 
attenuated by the application of a L-type Ca 2+ channel blocker |2j. 

In addition to the relative contribution of Ca 2+ influx to intracellular cal¬ 
cium levels, SCN neurons also have a rhythmically regulated reservoir of cal¬ 
cium that is not driven by membrane events. Both processes, inflow/outflow 
of calcium across the cell membrane and fixation/release of Ca 2+ from the 
calcium store, cooperate to generate intracellular Ca 2+ oscillations that seem 
to be crucial in driving a robust rhythm in gene expression. Indeed, rhythms 
in Ca 2+ levels, with its peak occurring during the day, seem to be a general 
feature of circadian systems [lj. 

It turns out that the oscillatory behavior of the intracellular concentration 
of Ca 2+ , and those of the variables of the genetic network as well, have a 
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much slower time scale than that of the variables involved in the generation 
of action potentials by the cell membrane. Thus, from a theoretical point 
of view, the problem is to understand how fast variables can control the 
dynamics of slower ones and vice versa. 

A conductance based model of spiking in cells of SCN has recently been 
presented in M- In those works, the huge variety of ionic currents that 
contribute to the membrane excitability was modeled in terms of the sodium, 
potassium and leakage currents characteristic of the Hodgkin-Huxley formal¬ 
ism plus a specific calcium current. They describe the periodic alternation 
between silent and excited states of the neuronal membrane in terms of direct 
and inverse Hopf bifurcations induced by the externally imposed time depen¬ 
dence of the ionic conductances. Nonetheless, in those works, no attempt is 
made to link this driving to changes in the intracellular level of either Ca 2+ 
or of any other substance. 

The aim of this paper is to suggest a plausible mechanism linking the 
electrical activity in the cell membrane to the circadian rhythms of the genes 
expressions, taking into account the two very different time scales associated 
to the membrane voltage (fast) and the genetic (slow) processes. The basic 
ingredient is to admit that the intracellular level of Ca 2+ controls both the 
firing frequency and the transcription of a clock gene. The control of the 
firing frequency is achieved by the free Ca 2+ concentration that, in turn, is 
controlled by the Ca 2+ in the cell cytoplasm coming from reservoirs. The 
reservoir Ca 2+ concentration varies in a much longer time scale than the in¬ 
tracellular Ca 2+ one. At the same time, the intracellular Ca 2+ concentration 
evolves slowly compared to the other variables characterizing the membrane 
electrical activity. On the other hand, the control on the circadian variables 
is established by assuming that the activation of the clock gene expression is 
proportional to the free Ca 2+ concentration. 

2. Ca 2+ controlled electrical activity 

We take the electrical activity taking place in the neuron membrane to 
be described by the two-dimensional system 


x — y — ax 3 + bx 2 + q — sz, 
y = c — dx 2 — y. 


( 1 ) 

( 2 ) 


The variable x represents the (dimensionless) voltage across the cell mem¬ 
brane and y is a recovering variable that describes the currents restoring the 
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polarity of the membrane after the emission of an action potential, q is a con¬ 
trol parameter. In the first of these equations, z stands for the free calcium 
concentration inside the cell, whose dynamics is supposed to be governed by 
the reactions 

Z — —$■ 0, 0 ) Z Zji ——$■ Z. (3J 

The first equation describes the outflow of Ca 2+ from the cell whereas the 
second one expresses the inflow of these same ions while the voltage-sensitive 
calcium channels are open during the membrane activation. Thus, the rate 
of this process is made dependent of the voltage. The last reaction describes 
the release of calcium ions from a intracellular reservoir described by the 
variable zr. The changing in Ca 2+ levels induced by the reservoir seems to 
be a critical part of the output pathway by which intracellular processes drive 
rhythms in neural activity [T]. The kinetic of calcium ions are completed by 
the reactions 

Z^Zr Zr^0 (4) 

where the first reaction describes the subtraction by the reservoir of free 
calcium ions and the second one corresponds to the direct loss of calcium 
ions by the reservoir. This set of reactions allows us to write the system of 
differential equations for the variables z and zr 


z = e(k(x) - k 2 z + k 3 z R ), (5) 

Zr = /x(k 4 z - k 5 z R ). (6) 

If we admit that the temporal scale for the evolution of the reservoir is 
much slower that the corresponding to the free calcium (// <C e), we have 

Zr = 0, =>- Zr = T, (7) 

T being constant on the e time scale. Thus, introducing this constant in Eq. 
[6] and assuming that k{x) = kix, we can write 

x — y — ax 3 + bx 2 + q — sz, (8) 

y = c-dx 2 -y, (9) 

z = e(k 1 x - k 2 z + g) , (10) 

with g = k 3 T. This three-dimensional dynamical system is formally identical 
to the Hindmarsh-Rose (HR) neuronal model [10]. Here, however, the slow 


4 


adaptation variable z introduced by these authors is reinterpreted as the 
intracellular free Ca 2+ concentration, whose level controls the much faster 
spiking variables x and y. The HR model shares essential qualitative features 
with several conductance-based models of excitable cells giving rise to square- 
wave bursting 13 Eld- All those models involve the interplay between a stable 
branch of stationary points and a limit cycle through a saddle middle branch 
in the equilibrium curve. By inducing the periodical transitions between 
these two coexisting attractors by means of a slow control variable, this 
dynamical mechanism allows the description of periodical silent and active 
phases in the neuron behavior. This same dynamical mechanism has recently 
been invoiced to interpret some experimental findings related to crustacean 
patterns generators [ 8 j. 

The HR system with e<l has two very different time scales, one of them 
associated with a fast subsystem (x and y variables) and the other associated 
with the dynamics of a slow subsystem (z variable). In fact, e 1 implies 
that in the time scale in which both x and y vary, z can be considered as 
a constant. Thus, calling sz = 7 we are led to analyze the dynamics of the 
two-dimensional system 

x — y — ax 3 + bx 2 + q — 7 , 

y = c-dx 2 -y , (11) 

in terms of the values of the parameter 7 . The equilibrium values x* of this 
two-dimensional system obey a cubic equation resulting from the crossing of 
both nullclines associated with Eq. (JTTj) , 

ax 3 + (d — b)x 2 — q — c + 7 = 0 . ( 12 ) 

Several regimes can be found as the parameter 7 is varied. As seen in Fig jT| 
for small negative values of 7 the only attractor is an equilibrium point (a 
stable focus) with non-zero x*-values. As 7 increases, the system undergoes a 
Hopf bifurcation rendering unstable the branch of equilibrium points giving 
rise to a stable limit cycle corresponding to voltage spikes (point B in FigjTJ. 
For still larger 7 values, the cubic equation has two unstable values coex¬ 
isting with a stable one. Furthermore, the limit cycle disappears through a 
homoclinic bifurcation when it collides with the unstable branch of x* (point 
A in the same figure). Notice the rather small interval of 7 values for which 
the stable limit cycle coexists with the stable equilibrium point appearing 
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through a saddle-node bifurcation. For still larger values of 7 , just a stable 
stationary branch remains. 

Following a technique pioneered by J. Rinzel some decades ago 112 a 
rather useful picture of the HR model dynamics as given by Eq. (10) can 
be obtained by projecting its orbits onto the xz-plane, and superimposing it 
with the bifurcation diagram for the fast subsystem described by Eq. (11). 
In Fig|2j an enlarged version of the coexistence region is presented including 
the nullclinc of the slow subsystem in ( 10 ) x = [k 2 z — g)/k\ (the red dot- 
dashed line). The solid black thin line in this figure corresponds to the 
bursting oscillations depicted in red in Fig. [T| This curve has been obtained 
from the x(t) and z(t) obtained by numerically integrating the full system 
of differential equations (10). It corresponds to the projection on the x — z 
plane of the global three dimensional attractor. 



Figure 1 : Z— shaped steady state curve x* = x*(y) for the fast subsystem (EqfTTj). 
The stable branch is depicted by a continuous line whereas the dashed lines repre¬ 
sent unstable branches of equilibria. The limit cycle appearing at 7 = —11.293 and 
ending at the homoclinic bifurcation at A appears in gray. The projection of the 
global attractor of the full three-variable system (in red) has been superimposed to 
this diagram. Parameter values are a = 1.0, b = 3.0, c = 1.0, d = 5.0, q = 0.3, k\ = 
1.0, k 2 = 0.1 and g = 1.33. 

As one can observe in Fig. [ 2 j when z slowly moves to the left (notice 
that z < 0 under the z—nullcline), the x variable evolves along the branch 
of stable equilibria (C -A D ) until it comes to an end at a saddle-node 
bifurcation (point D ). At this point it is forced to jump to the limit cycle 
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Figure 2: Enlarged version of coexistence region depicted in Fig. [l] The homo¬ 
clinic bifurcation occurs at A and the fold bifurcation at D. The red dot-dashed 
line corresponds to the z—nullcline. The solid black thin line corresponds to the 
bursting oscillations depicted in red in Fig. [1} 


by crossing the z—nullcline and afterwards it is forced to move to the right 
(now z > 0) while performing fast oscillations. This oscillatory behavior ends 
when the limit cycle disappears through the homoclinic bifurcation at A. At 
this point, the x variable must jump to the stable branch of equilibria and 
the whole process repeats itself. 

The temporal evolution of x(t) and z(t) has been depicted in Figj3j As 
we can see, the free running x(t) variable alternates periodic episodes of fast 
spiking with silent epochs of much slower evolution. On the other hand, 
z(t) evolves in the slower time scale. In this case, the temporal scale of the 
transitions between the silent and the spiking states of the neuron coincides 
with that of the z variable. Each period of the z oscillations corresponds to 
a cycle characterized by a rise of the calcium level associated with the burst 
of spikes and followed by a silent phase in which the concentration of Ca 2+ 
slowly decays. 

Note that our model describes the membrane activity without explicitly 
incorporating the ionic currents. Nonetheless, our equations show a bifurca¬ 
tion structure close to that of more complicated conductance models [3]. 
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3. Coupling electrical and molecular activities 

On the genetic side of the model, we have considered the transcription 
of only one clock gene. Thus the variables of interest are the intracellular 
concentration of clock gene mRNA, the corresponding level of a clock protein 
and that of a transcriptional inhibitor whose action closes the self-regulatory 
feedback loop. All these variables obey the differential equations of the well- 
known Goodwin model [9j. 

To link the membrane electrical activity to the clock gene expression, we 
assume that the level of free calcium acts as an activator of the transcription 
of the mRNA. Then, our complete model is embodied in the six-variable 
system of differential equations 


X = 

y — ax 3 + bx 2 — sz + q + pY, 

(13) 

y = 

c — dx 2 — y, 

(14) 

z = 

e(k 1 x - k 2 z + g), 

(15) 

X = 

e( aZ kx), 

(16) 

Y = 

e(k f X — kY), 

(17) 

Z = 

e(k f Y - kZ) 

(18) 


where A", Y and Z describe the clock mRNA, the clock protein and the 
inhibitor of the transcription, respectively. Throughout this work we will 
take a = 1,6 = 3, c = l,d = 5,s = 1, h — 10, k = 2, kf = 2 and, a = 8. 
Note that the above equations contemplate a direct feedback mechanism by 
which the molecular clock is able to drive the electrical activity of the cell 
membrane through the term pY appearing in the voltage equation. At the 
same time, the variable 0 representing the free Ca 2+ concentration inside the 
cell, influences the clock gene concentration X through a Hill-like term. For 
more elaborated circadian models involving several clock genes [131 E] the 
coupling between the membrane and genetic activities could be formulated 
along the lines of our methodology by including terms describing the Ca 2+ 
activation of the different gene expressions. 

We first study the case with no feedback from the genetic variables on 
the membrane ones (p = 0). In Fig. [3] we have depicted the temporal 
behavior of the voltage variable and that of the intracellular concentration 
of free calcium as well as the evolution of all the variables of the genetic 
subsystem. The parameter e sets the spiking frequency of the system so that 
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Figure 3: Square-wave bursting of the model with p = 0 as a function of (dimen¬ 
sionless) time et (upper panel).The evolution of X(t) (black), Y(t) (blue) and Z(t) 
(red) as well as the temporal behavior of z(t ) (red thick line) have been depicted in 
the lower panel. Parameter values are e = 0.001, p = 0, q = 0.3, k\ = 1.0, k ‘2 = 0.8 
and g = 1.23. 


the smaller it is, the higher the frequency of spiking achieved during the 
bursting phase. The value of q is critical to the character and persistence of 
these bursts. In order to save computation time we have set the values of 
e small but with no intention to fit quantitatively the experimental values 
for the neurons of the SCN. One could interpret the behavior observed in 
Fig. | by noting that the start of a burst of spikes activates the inflow of 
Ca 2+ , thus increasing its intracellular level until it forces the membrane back 
to its resting electrical state. On the other hand, the slow oscillation in the 
intracellular calcium concentration forces a slightly delayed oscillation in the 
level of mRNA (black line in the lower panel), leading to an additional delay 
in the evolution of the protein level (blue hue) as well as in the inhibitor (red 
line). Note that there is no modulation of the bursts amplitude due to the 
fact that we have chosen p — 0. 

Now we turn to the case with p ^ 0. As we can see in Figj4j, the existence 
of the pY feedback term in the first equation of the model generates a rel¬ 
atively complex temporal modulation of the voltage variable along the day. 
Again, the genetic variables oscillate at the slower time scale of the voltage 
variable and with a certain dephasing among them. As we can see in Fig|4j 
the level of the clock gene mRNA peaks at the start of the silent phase of 


9 
































the neuron that takes place during the daylight hours (see Figj5| while it 
decreases during the night. 



Figure 4: Autonomous behavior of x(t) and z(t) for the model with p / 0 (upper 
panel). The evolution of X(t) (black), Y(t) (blue) and Z(t ) (red) have been 
depicted in the lower panel. Parameter values are k\ = 4, A ^2 = 0.6, q = 0.5 and 
p = 12.5. 

In Fig|5]we present the autonomous behavior of the voltage variable in 
terms of the Zeitgeber Time (ZT) during a whole day. The time origin for 
the data plotted in this figure has been shifted to coincide with that of the 
data reported by Belle and coworkers for the forcing of perl-containing SNC 
neurons with a sequence of equal periods of light and darkness [Tj. The results 
presented here are qualitatively close to those experimental findings although 
some differences in timing remains. From ZTO to ZT2 approximately, the 
neuron is firing, whereas from ZT4 to ZT10, it has been deeply depolarized 
and has become silent. The value chosen for p in the equations of the model 
affects the period of depolarization. For p = 13.5, the spiking behavior starts 
again at ZT10 and lasts until ZT18, when the neuron becomes hyperpolarized 
and its activity ceases. This silent period increases as the parameter 
decreases. Finally, at Z22 the neuron starts firing again. The intracellular 
concentration of calcium ions, on the other hand, rises during the day, peaks 
in the evening and early night and decreases to a minimum along the late 
night hours. Note that the spiking frequency changes with the ZT value. 
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Figure 5: Autonomous behavior of the voltage variable (upper panel) and that of 
the concentration of intracellular calcium (lower panel) during a whole day. The 
silent phase characterized by a strongly depolarized membrane appears through 
an inverse Hopf bifurcation and ends at a direct bifurcation of the same kind. 
Parameter’s values are k\ = 4, = 0.6, q = 0.5, and p = 13.5. The time origin 

has been shifted to coincide with that used in [J]. 


4. Concluding remarks 

Experimental evidence indicates that the link between electrical activity 
and clock gene expression in SCN cells is provided, at least in part, by the in¬ 
tracellular dynamics of Ca 2+ ions associated with the existence of a reservoir 
of these ions inside the cell. The nature and localization of this intracellular 
reservoir is unknown at present although it seems that this role is played by 
the endoplasmic and sarcoplasmic reticulae [ 15] . 

In this work we have presented a mathematical model aiming to explore 
a simple dynamical mechanism for the driving of clock genes by the firing 
activity of neurons. Some work carried out in the last few years suggests the 
existence of a control of the bring frequency by clock gene expression, but 
very little is known about the mechanisms by which genetic oscillations are 
able to drive rhythms in neural activity [TJ. So, in this work we have explored 
this issue by assuming that the protein produced by the transcription of the 
clock gene modulates the neuronal bring. Similar results would have been 
obtained by using the mRNA variable or the inhibitor to control the neuronal 
bring. The role played by a calcium reservoir has been taken into account by 
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considering the kinetics of calcium fixation and its release by the reservoir. 
It seems possible to extend the main ideas that we have elaborated here to 
link the genetic and membrane variables of the SNC using a conductance 
based model for the membrane activity as well as a dynamics of more than 
one clock gene. Work along these lines is in progress. 
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